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In systems belonging to the universality class of the random field Ising model, the standard hyper- 
scaling relation between critical exponents does not hold, but is replaced by a modified hyperscaling 
relation. As a result, standard formulations of finite size scaling near critical points break down. 
In this work, the consequences of modified hyperscaling are analyzed in detail. The most striking 
outcome is that the free energy cost AF of interface formation at the critical point is no longer a 
universal constant, but instead increases as a power law with system size, AF oc L° , with the 
violation of hyperscaling critical exponent, and L the linear extension of the system. This modified 
behavior facilitates a number of new numerical approaches that can be used to locate critical points 
in random field systems from finite size simulation data. We test and confirm the new approaches 
on two random field systems in three dimensions, namely the random field Ising model, and the 
demixing transition in the Widom-Rowlinson fluid with quenched obstacles. 



I. INTRODUCTION 

Understanding the effects of quenched random disor- 
der onphase transitions has been a longstanding chal- 
lenge 16]. Analysis of experiments on such systems is 
typically more difficult than work on pure systems Q. 
Theoretical methods are hampered by the fact that, for 
spin glasses and systems exposed to random fields, the 
marginal dimension d* = 6 (the marginal dimension is 
the dimension above which mean field theory is believed 
to be reliable). In contrast, for pure systems, d* = 4 
0,0- As a consequence, predictions of renormalization 
group expansions in e = d* — d dimensions tend to be less 
reliable in the physically relevant dimensions (d = 2,3) 
when quenched disorder comes into play. Computer sim- 
ulations, albeit very useful for the study of critical phe- 
nomena in pure systems jjj-llljl. suffer from the problem 
that for systems exhibiting quenched random disorder 
an additional average over many samples drawn from 
the distribution characterizing the disorder needs to be 
taken. The disorder average, denoted [•], comes in addi- 
tion to the usual thermal average, denoted (•), and hence 
the computational effort is order of magnitudes larger. 
Since most analysis of critical phenomena by simulations 
[9l-fli1| relies on finite size scalin g ( l2l-[l6j. lack of self- 
averaging in random systems [17H20| is also a problem. 

For Ising ferromagnets diluted with nonmagnetic im- 
purities there is no doubt that the transition, from the 
high-temperature disordered to the low-temperature or- 
dered phase, remains second order in d = 2, 3 Q. In 
addition, the hyperscaling relation Q between critical 
exponents remains valid 



2-a = 2/3+j = du (d = 2,3), 



(1) 



and rather accurate estimates for these exponents are 
available 2l| (we use standard symbols to denote the ex- 



ponents; definitions arc provided in Scction[n]). For Ising 
ferromagnets in random fields, however, the situation is 
radically different. In d — 3 dimensions, the existence of a 
transition at nonzero temperature was controversial until 
a proof for the existence of a spontaneous magnetization 
settled this issue [22j ; rigorous results on the order of the 
phase transition in the d = 3 random field Ising model 
(RFIM) are still lacking however. If one accepts the evi- 
dence from numerical studies [HI, HiJ that the transition 
is second order, it must have very unconventional criti- 
cal behavior [25l - l27] . The key point is that the standard 
hyperscaling relation, Eq. (JlJ, no longer holds, but is re- 
placed by 



2 - a = 2/3 + 7 = v(d - 



(2) 



The exponent 9, called the "violation of hyperscaling ex- 
ponent" , is believed to be [28M30| 



j/u = 2 — 7], 



(3) 



where the critical exponent r\ describes the decay of the 
spin pair correlation function right at the critical temper- 
ature 0. In addition, it is believed [25l427| that critical 
slowing down in the RFIM is not described by the usual 
power law for the relaxation time r a with z the 
"dynamic critical exponent" [3lT | and £ the correlation 
length, but instead is governed by a much more severe 
"thermally activated critical slowing down" [25U27} 



lnr cx £ fl 



(4) 



with T — > T c from above. It should not come as a sur- 
prise that Eqs. ([2]), © and ((3]) make the study of the 
RFIM by Monte Carlo (MC) methods very difficult, and 
early studies even claimed a weak first order transition 
[32l | . Another problem is that standard finite size scaling 
formulations typically rely on the validity of hyperscal- 
ing, which does not hold in the RFIM [§, Q. Some of 
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the consequences resulting from the violation of hyper- 
scaling, Eq. @, were already noted in previous works 
[33j . and exploited in recent studies of colloid-polymer 
demixing in random porous media @. 

The purpose of the present paper is to analyze the con- 
sequences of Eq. ([2]) for finite size scaling in more detail. 
In particular, we shall focus on the free energy barrier 
AFl separating the coexisting phases for T < T c in a fi- 
nite system of linear extension L. It has been found that 
AF^ increases quite strongly with L at T = T c in the 
RFIM [IB)]. This behavior is puzzling because a grow- 
ing barrier is usually associated with a first-order phase 
transition. In this work, the theoretical justification for 
this behavior is provided. We will show that the barrier, 
which in the regime where the transition is first-order 
scales as 



AF L = 2/ int L^ i (T<T C ), 



(5) 



with /i n t the interfacial tension |15| . right at the criti- 
cal point is related to the hyperscaling violation critical 
exponent 



AF L ocL e (T = T C ). 



(6) 



The factor-of-two in Eq. is a consequence of periodic 
boundary conditions, which induce two interfaces in the 
system when T < T c . We will provide numerical evidence 
in favor of Eq. ([6]) using simulation results obtained for 
two random field systems in d = 3 dimensions, namely 
the RFIM and the Widom-Rowlinson fluid with quenched 
obstacles. We emphasize that in d = 2 dimensions the 
RFIM is without a phase transition, in which case the 
analysis of the present paper does not apply. 



II. THEORETICAL BACKGROUND 

We consider a system of N Ising spins, situated on a 
d-dimcnsional lattice of linear extension L with periodic 
boundaries, inside an external magnetic field H. The 
instantaneous magnetization per spin is defined as 



1 N 

5?> 



N 



(7) 



with Si = ±1 the value of the spin at the z-th lattice site. 
We assume that the system, in the thermodynamic limit 
L — > 00, undergoes a second-order phase transition at 
critical temperature T = T c and field H = H c . Following 
standard practice, we introduce the relative deviations 



t = T/T c -l, h = H/H c -l. 



(8) 



In the vicinity of the critical point (t = 0,h = 0), the 
specific heat C, susceptibility x, and correlation length £ 
diverge as power laws 



C(t,h = 0) oc \t\- a , x(t, h = 0) oc |i| 
£(t,fc = 0) cx \t\-". 



(9) 



In the ordered phase T < T c , a finite magnetization M 
(order parameter) and interfacial tension f- lnt develop 



M(t,h 
fint(t,h 



0) cx \tf 
0) cx |*|" 



(t < 0), 
(t < 0). 



(10) 



We first give a heuristic derivation of the standard hy- 
perscaling relation, Eq. ([T]), which is valid in pure sys- 
tems, i.e. without quenched random fields. Following the 
static scaling hypothesis (H) , the singular part of the free 
energy density takes the form 



Mt,h) = \t\ i -°f(h/\tf+r), 



(11) 



with f(x) a scaling function. The order parameter is 
obtained by differentiating / s i n g once with respect to the 
external field 



M(t,h = 0) oc 



dh 



« 1*1 



2— a—/3— 7 



(12) 



h=0 



which, upon comparing to Eq. (|10[) , immediately yields 
2 — a = 2/3 + 7 Hi- Near criticality, the singular part of 
the free energy can be attributed to correlated regions of 
spins (clusters) of linear dimension £ 0, H3] • Each cluster 
has essentially one Ising degree of freedom (magnetiza- 
tion direction up or down), and can orient independently 
from its neighbors. Thus, while at T — > 00 the total 
free energy F of the system is due to the entropy of N 
non-interacting spins, F = — (fcsTln2)A^, near T c we can 
attribute the singular part of F to the entropy of N/£ d 
independent clusters of spins, F = — (fc^T In 2)iV/£ d , and 
hence 



/ sing M)c<r d oc|*r 



(13) 



where in the last step Eq. ([9]) was used. Comparing the 
above equation to Eq. (fTTj). the standard hyperscaling 
relation 2 — a = dv immediately follows. 

For an Ising system in quenched random fields near 
criticality the situation is different. To be specific, con- 
sider a random field ±r acting on each spin (with the 
signs ± drawn with equal probability such that [r] =0). 
We can still split the system into clusters of linear dimen- 
sion £, such that each cluster may be considered as inde- 
pendent from its neighbors. However, the main contribu- 
tion to /sing in this case is not the entropy, Eq. (|13[) . but 
rather the Zeeman energy due to the coupling to the ran- 
dom field. In a region of volume £ d , the sum of the ran- 
dom fields exhibits Poissonian fluctuations ±r£ d / 2 . The 
random field excess per spin is therefore of order 



Ar oc ±r£ 



-d/2 



(14) 



which may be conceived as an external field acting on the 
spins in the region. This implies a finite magnetization 
per spin 



(m) oc X Ar oc ±r X £, d/2 , 



(15) 
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with x the susceptibility. The Zeeman contribution to 
the free energy thus becomes 

f sing (t, 0) cx (m)A R cx r 2 X ^ d <* l^" 7 , (16) 

which dominates the entropy contribution, Eq. (|13[) . 
upon approach of the critical point t — > 0. If we insist 
that / s i n g retains the scaling form of Eq. (fTTj) , it follows 
that 2 — a = dv — 7; using 6 = 7/1/ (Eq. ((3])) then yields 
the modified hyperscaling relation (Eq. ([2])). 

Next, we consider the exponent /i of the interfacial 
tension (Eq. (|10[0 . Since Eq. (|TT|) is a free energy per 
volume, and since near T c the correlation length is the 
only relevant length scale, a simple dimensional argument 
implies that 



/int °C /sing£ ~>fJ> = 2- a-V. 



(17) 



In case of hyperscaling this implies /Lt = (d — \)v. The 
important point of the present discussion is that this re- 
lation does not hold in the RFIM, since the hyperscaling 
relation, Eq. (TTJ), is violated and replaced by the mod- 
ified relation, Eq. ©. We can still infer that Eq. ([17]) 
should hold, but now one must use Eq. @, which leads 
to /xrpim = (d—1 — 6)v. Finally, we discuss fluctuations, 
which are typically large near phase transitions. In the 
pure Ising model, the thermally averaged magnetization 
plays the role of order parameter M, while the suscepti- 
bility x reflects its thermal fluctuations 



pure Ising model 



M = (|m|), 



X = L*{(m 2 )-(\m\) 2 ) 
For the RFIM, the obvious generalizations are 



RFIM 



M=[(\m\)}, 

Xcon = L d [(to 2 ) - (H) 5 



(18) 



(19) 



with [•] the disorder average (factors of ksT have been 
dropped in our definitions), x an d Xcon are called "con- 
nected" susceptibilities: they reflect thermal fluctua- 
tions, which are present in both models, and diverge at 
criticality with exponent 7 (Eq. ((9])). Note that our defi- 
nitions of the order parameter and susceptibilities use the 
absolute value of the instantaneous magnetization, as is 
commonly done in simulations (36j ] . 

In the RFIM, we can also define a "disconnected" sus- 
ceptibility 



27] 



RFIM xdis = L d ([{\m\)< 



(20) 



which docs not have its analogue in the pure model (re- 
moving the disorder average [•] trivially yields Xdis = 0). 
The motivation to introduce Xdis stems from the obser- 
vation that {\m\) depends on the random field sample. 
Hence, in the disorder average, there will be sample-to- 
sample fluctuations in (|m|), which is precisely what Xdis 
corresponds to. Upon approach of the critical point, the 
disconnected susceptibility also diverges 



TABLE I. Critical exponents of the 
RFIM in d = 3 dimensions; see Ref. 
list of results for the RFIM. 



ure Ising model and 
for a more elaborate 



pure Ising 



RFIM 



0.326 [8] 
0.630 [8] 

1.240 [8] 



0.0-0.02 [37J 
0.06 [24J 
[23] 

1.14 [37], 1.67 [37] 

1.02 [24] 

1.1 [23] 

2.25 [38] 

1.9 [24] 

3.4-5.0 [37] 

2.9 [24] 



defining a new critical exponent 7. It is predicted that 
7 = 27 [25l - [30l ] . implying that sample-to-sample fluctu- 
ations dominate over thermal ones at criticality. If we 
substitute, in Eq. ©, 6 — > j/i/ and 7 — > 7/2, the modi- 
fied hyperscaling relation becomes 



2/3 + j = dv, 



(22) 



which is just the standard hyperscaling relation, Eq. ((T|), 
but with 7 replaced by 7. 



III. FINITE SIZE SCALING 
A. pure Ising model 

We first consider finite-size scaling (FSS) in the pure 
Ising model in d = 3 dimensions. The Hamiltonian is 
given by 



H 



Ising 



= -jJ2SiSj-Hj2Si, J>0, 



(23) 



Xdis oc |i| 7 , 



(21) 



with a sum over nearest neighbors. In what follows, 
the temperature T is expressed in units of fcs/J, with 
ks the Boltzmann constant. For the d = 3 Ising model 
on cubic periodic lattices, the critical temperature T c sa 
4.511 [HI. The critical exponents are known relatively 
precisely, although not exactly (Table [XJ> . 

A key quantity in the numerical study of phase tran- 
sitions is the order parameter distribution (OPD), de- 
noted Pi (to), and defined as the probability to observe 
the system in a state with magnetization m. We assume 
the OPD is normalized: P L (m) dm = 1. The OPD 
depends on the system size L, and on the control pa- 
rameters T and H. In the pure Ising model, due to spin 
reversal symmetry, the critical field H c = 0, and so we 
set the external field to zero. The OPD is then an even 
function, Pl{— m) = Pl(to), irrespective of T and L. 

In the ordered phase, T < T Cl the transition is first- 
order. There exists a spontaneous magnetization, which 
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may be positive or negative. The OPD is a superposition 
of two Gaussians, centered around to = ±m , with expo- 
nentially small finite size effects in the peak positions [39| ■ 
The definition M = (|m|) corresponds to (half) the peak- 
to-peak distance. In the disordered phase, T > T c , the 
OPD tends to a single Gaussian peak centered around 
m = 0. In both cases, the system self-averages: the 
peak widths decay oc L~ d / 2 , ultimately becoming sharp 
(5-functions. 

At criticality, the L-dependence of the OPD is given 
by the scaling form 



(a) 1.0 



-L-oepf 

Mm 



P L (m) oc P{L"m) -> (m k ) oc iT fc6 (T = T c ) (24) 



with b a constant, and P(x) a scaling function charac- 
teristic of the Ising universality class. The standard FSS 
expressions for the order parameter and susceptibility are 

EMI 



M oc L- p /'\ X oc L 7/ " (T = T c ). 



(25) 



In order to be consistent with these expressions, Eq. (f2~4")) 
requires b = (3 /is, and the validity of standard hyperscal- 
ing. For the pure Ising model, the OPD at criticality is 
bimodal with overlapping peaks (Fig. UKa)). If one plots 
the distributions versus the scaling variable, x = L^' v m, 
the curves for different L overlap (Fig. [TJb)). A further 
consequence of Eq. ([24]) is that cumulant ratios such as 



2\2 



XJ X = (m*)/(\m\y, Ui = (m 4 )/(m 2 ) 



etc. 



(26) 



are L-independent at criticality. Since the peaks in the 
critical OPD of the pure Ising model overlap (as opposed 
to being sharp), the corresponding cumulant ratios are 
distinctly different from the off-critical values. For ex- 
ample, considering the U\ cumulant, it holds that 



lim Ui = 




(27) 



where Uf = 1.2391(14) for the d = 3 Ising model [41 
This behavior is extremely useful to extract T c from sim- 
ulation data [9I4TTI H3 [. see Fig. [Ijc). Note that Ui is 
essentially the ratio between the order parameter and its 
thermal fluctuations 



2 _ 



(to 2 ) — (|m|) 5 

<H> 2 



= Ux-l. 



(28) 



The fact that U\ = Iff at criticality implies that ctt re- 
mains finite. Put differently: the thermal fluctuations 
of the order parameter M remain comparable to M it- 
self at T c . From this consideration one also understands 
why the peaks in the OPD are broad and overlapping. 
Alternatively, we may write 



(b) 1.0 



(c) 1 .45 




4.58 



FIG. 1. FSS in the d = 3 pure Ising model in the critical 
regime. The linear dimension L is given in units of the lattice 
spacing, (a) OPD P L {m) obtained at T = T c , H = H c = 0, 
and for several system sizes, (b) The same data plotted versus 
the scaling variable, with the critical exponents taken from 
Table UI the data for different L overlap, (c) Demonstration 
of the cumulant intersection method to locate T c . Plotted is 
Ui versus T for several system sizes. At the critical point, the 
curves for different L intersect. In the ordered (disordered) 
phase, U\ —¥ 1 (U\ —¥ tt/2) as L increases, in accord with 
Eq. (27]). 



4 = X /L d M\ 



(29) 



-1.0 -0.5 0.0 0.5 1.0 

m 

FIG. 2. Finite size effects in the d = 3 pure Ising model in the 
ordered phase, where the transition is first-order. Plotted is 
the scaled-and-shifted logarithm of the OPD, Wl = In Pi, (m), 
for various system sizes at T = 3.33, which is well below T c . 
The peak height corresponds to the interfacial tension /i n t. 
Note also the flat region unfolding between the peaks as L 
increases. 



from which, using Eq. (|25[) and hyperscaling, one also im- 
mediately derives that <tt is L-indepcndcnt at criticality. 

The free energy barrier is obtained from the logarithm 
of the OPD W L = In P L (m). We define AF L as the 
average peak height, measured from the minimum "in- 
between" the peaks. In the ordered phase, T < T c , the 
transition is first-order, and the barrier is related to the 
interfacial tension /; nt via Eq. ((SJ). This is shown in Fig. [5] 
Note also that the peak positions - at least on the scale 
of the graph - do not reveal any strong L dependence 
either, consistent with exponentially small finite size ef- 
fects [39j . The peaks also become sharper as L increases, 
showing that the system is self-averaging. Finally, we 
note that a flat region between the peaks in Wl unfolds 
as L increases. This is a sign that interactions between 
the interfaces through the periodic boundaries are van- 
ishing H3. 

Precisely at criticality, the scaling of the barrier is dif- 
ferent. We may still assume that Eq. ((5j) holds, but on a 
length scale that is set by the correlation length 

AFjt oc /irtf 1 - 1 oc ^-Ma-"-")/", (30) 

where in the last step the critical power law of /i n t and 
Eq. (|17[) were used; by virtue of hyperscaling, the length 
scale drops out. Hence, the barrier is a constant L- 
independent value AF^ = AF* at criticality. Of course, 
the fact that the OPD at T = T c has a universal shape, 
see Fig. QJb) , also implies this property. In the pure Ising 
model, the barrier thus scales as 

[2/i nt L d - 1 F<T C , 
lim AFl = < AF* T = T c , (31) 



4.44 4.46 4.48 4.50 4.52 4.54 4.56 4.58 
T 

FIG. 3. Scaling of the free energy barrier AFl in the pure 
Ising model. By plotting AFi versus T for various L, the 
critical temperature appears as an intersection point. The 
horizontal line marks the (universal) value AF* ~ 0.9 for the 
d = 3 Ising model. 



This behavior is also well-suited to locate T c [43J. For 
instance, one plots AFl versus T for various L; at the 
critical point, the data for different system sizes intersect 
(Fig.EJ). 

In brief, we have summarized FSS in the pure Ising 
model. The important message is that, due to hyperscal- 
ing, the OPD assumes a universal shape at the critical 
point. As a result, the free energy barrier and selected 
cumulant ratios assume non-trivial L-independent values, 
which can be used to locate the critical point. We also 
note that, by using the intersection methods of U\ and 
AFl (Fig.QJc) and Fig. [3]), moderate system sizes suffice 
to locate the critical point with an accuracy better than 
one part in a thousand. 

B. RFIM 

We now consider FSS in the RFIM at its critical point. 
The Hamiltonian reads as 

Hrfim = -J S i S J - ^ r iSi - H^Si, (32) 

<i,j> i i 

J > 0, with Ti the quenched random field acting on the 
spin at the i-th lattice site. It is convenient to draw Vi 
from a distribution that is symmetric about zero; this 
ensures that H c = 0. The most common choices are the 
bimodal distribution F(r^) oc 8(ri — a) + S(ri + a) and the 
Gaussian F(r^) oc exp(— r|/2tr 2 ), where a is the random 
field strength. In contrast to the pure Ising model, the 
critical exponents of the RFIM are not known very pre- 
cisely, see Table HI where several exponent estimates from 
theoretical and simulational works are listed. We believe 
it is safe to conclude that (3 is close to zero. Modified 
hyperscaling then implies 9 = 7/1/ « 1.5 and j/v ~ 3, 
where dimensionality d — 3 is assumed. 
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1. Consequences of modified hyperscalmg 



For instance: 



One of the most striking consequences of hyperscaling 
violation is that the thermal fluctuations become negli- 
gible at the critical point. For the RFIM, the analogue 
of Eq. (HU becomes a\ = x con / L d M 2 ; using the FSS 
expressions M oc L~^l v and Xcon oc L^/ u , it follows that 



o£ oc L 2 ^ v+ "i' v - d oc L-^' v 



(T = T c ), (33) 



where now the modified hyperscaling relation, Eq. ©, 
was used. In the RFIM, ctt thus decays to zero with in- 
creasing L, whereas in the pure Ising model ctt saturates 
at a finite L-independent value. 

Even though the thermal fluctuations vanish in the 
RFIM for large L, we must not forget about the samplc- 
to-sample fluctuations, which are characterized by Xdis- 
In line with ctt, we compare the order parameter to the 
magnitude of sample-to-sample fluctuations as 



Xdi: 



L d M 2 



oc L 2fi /' y+ ^"- d (T = T c ), (34) 



where also the FSS expression Xdis oc lHl v was used. 
The remarkable consequence of modified hyperscaling, 
Eq. is therefore that ctd oc L°, i.e. becoming con- 

stant at criticality. Hence, in the RFIM, it is the sample- 
to-sample fluctuations that "scale with L" , not the ther- 
mal fluctuations. 

How does this modified scaling affect the OPD? First 
note that, in addition to T, L, and LT, the probability 
to observe a certain instantaneous magnetization m also 
depends on the random field sample. We therefore write 
Pi i i(m), where the index i denotes one particular sam- 
ple of random fields. Wc thus have a set of distributions. 
In practice, this requires that Pl^to) be measured for 
at least i = 1,...,K samples, where K must be large 
enough. We can immediately rule out that Pi^m) at 
criticality obeys the scaling form, Eq. (|24|). since hyper- 
scaling is violated. Assuming that the majority of distri- 
butions Pi,,i(jn) remains bimodal at T c - which needs to 
be verified in practice - the peak-to-peak distance scales 
as the order parameter M, while the squared peak widths 
W 2 oc Xcon/L d - Since ot(= W/M) decays to zero, see 
Eq. (|3"3")l , it follows that the peaks in PL,i(m) become 
sharp. Again, this is in contrast to the pure Ising model, 
where the critical OPD features broad and overlapping 
peaks. In the RFIM, the shape of PL,i{m) at T = T c and 
T < T c is the same: bimodal with sharp non-overlapping 
peaks. The crucial difference is that, at T = T c , the 
peak-to-peak distance decays oc L~^l v ', while for T <T C 
the peak positions saturate at finite values ±mo- In 
the disordered phase, T > T c , PL,i{rn) should again 
be single peaked. As a consequence, ratios of connected 
quenched-averaged moments, such as [{m 2k )]/[(m k ) 2 } or 
[(m 2k )]/[(m k )} 2 , no longer assume "special" values at 
criticality, but equal those of the ordered phase T < T c . 



lim U - [(m2)] 

iSlo 1,con = KHp 




(35) 



which does not lend itself well to extract T c from finite- 
size simulation data. 

Does this imply there is no "scaling" at all in the 
RFIM at its critical point? The answer to this ques- 
tion is an unequivocal "No" ! Scale invariant distributions 
and observables still exist in the RFIM, but they must 
be constructed keeping modified hyperscaling, Eq. ([22]) . 
in mind. For instance, to each random field sample i 
there corresponds a distribution Pj Jj j(m), from which an 
average magnetization (|m|)j can be obtained. Due to 
sample-to-sample fluctuations, the values (|m|)j will gen- 
erally differ. Hence, it is useful to consider the distribu- 
tion ^((ItoD), defined as the probability of a particu- 
lar random field sample yielding a thermally averaged 
magnetization In the absence of quenched disor- 

der, "P_l((|to|)) reduces to a 5-function; in the presence of 
quenched disorder, Pl((|?tj|}) may retain a finite width. 
The moments of correspond to [(|m|) fc ], which 

are precisely the quantities needed to compute the or- 
der parameter and the disconnected susceptibility. If we 
compare the average of T > L({\m\}) to its root-mean-squarc 
width, we recover od of Eq. (J34J) ; by virtue of modified 
hyperscaling, the latter becomes constant at criticality. 
Hence, in the RFIM, it is the distribution VL((\ m \)) that 
remains broad at criticality. Our "Ansatz" is therefore 
that the scaling of the OPD in the pure Ising model, is 
replaced by scaling of 7- > l((|«'i|)) hr the RFIM. Wc thus 
propose 



r L ({\m\)) ^V(L\\m\)) 

[{\m\) k \ oc L- kb (T = T c , RFIM) 



(36) 



as the analogue of Eq. (|24p . with V(x) a scaling function 
characteristic of the RFIM. Consistency of Eq. (|36t with 
M oc L~ fi / y and Xdis oc IT*! 11 requires b = f3/v, and the 
validity of modified hyperscaling, Eq. ([22)1 . Note that 
Eq. ([5d| also implies that "disconnected" cumulants such 
as 



l(\m\) 2 }/[{\m\)f 



(37) 



become L-independent at T c . This property suggests 
that a generalization of the "cumulant intersection 
method", Fig. HJc), is also feasible in the RFIM. In this 
case, one should plot Lq.di s versus T; curves for different 
L should intersect at T c . 

The second main consequence of modified hyperscal- 
ing concerns the scaling of the free energy barrier. The 
barrier is no longer constant at T c , but rather AF^ oc L e , 
with 9 = 7/1/ the "violation of hyperscaling" exponent. 
This follows immediately from Eq. (|30[) , where now the 
modified hyperscaling relation must be used, as well as 
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(0,0) x 



FIG. 4. Schematic plot of the scaling function F(x), defined 
by Eq. (|39|l. describing the free energy barrier AFl in the 
RFIM versus the scaling variable x = tL 1 ^" . There occurs 
a smooth crossover from F(x) oc |^ | ^ d ~ i — e) ^ f Qr x <g. 0, to 
F(x) cc l/\x\ p for x > (with p > vff). 

the FSS "Ansatz" £ oc L [l2|-[l4|. We thus expect, for 
random field Ising universality, 

(fintL*- 1 T<T C , 
lim AFl oc { L e T = T c , (38) 

[0 T>T C . 

Following standard FSS practice fl2rll4l |. we may also 
write 

AF L = L 9 F(x), x = tL 1 ^, (39) 

with F(x) a scaling function. The scaling of the bar- 
rier in the ordered phase, T < T c , implies that F(x) oc 
^(d-i-e)v ^ x ^ Precisely at criticality x = 0, 

we should recover Eq. (j6]), i.e. F(0) > 0, while in 
the disordered phase T > T c the barrier should vanish 
F(x) cc l/\x\ p , p > v6 (i» 0). From these consider- 
ations, as well as from the fact that the scaling function 
must be smooth, we derive the sketch shown in Fig. U) 
The fact that F(x) is a smooth function, implies that 
(huge) free energy barriers AFl oc £ e persist above T c 
also (in sharp contrast to the pure Ising model) . The lat- 
ter give rise to the Arrhenius law for the relaxation time, 
Eq. g]). 

2. Practical considerations: tuning the external field 

In FSS studies, the critical region is "scanned" by vary- 
ing the control parameters T and H. Mathematically, 
this can be conceived as following a path in the (T, H ) 
plane. One may choose the path freely, as long as it 
passes through the critical point (T c , H c ) in the thermo- 
dynamic limit. In the pure Ising model and RFIM, the 
critical field H c = 0, and so the critical region may be 
scanned by varying T at fixed H = 0. We call this the 
symmetry path Is : H = 0. It may also happen that H c 



is not known beforehand. This is often the case in flu- 
ids, where the analogue of H is the chemical potential. 
In these situations, different paths must be constructed. 
One example is the "equal- weight" path [44| , whereby H 
is tuned such that the peaks in the OPD have equal area. 
The field now becomes a non-trivial function of temper- 
ature and system size H = f(T,L). As it turns out, an 
infinite number of paths can be constructed along these 
lines |45j . Here, we will mostly use the path lr, whereby 
H is tuned such that 

l r : d(m)/dH -> max. (40) 

Note that, when lr is used in conjunction with quenched 
disorder, H not only depends on T and L, but also 
on the random field sample i, that is H = fi(T,L) 
with i = 1,...,K. For fixed T and L, each sample 
thus yields its own field Hi. A sharp transition re- 
quires that, for T < T c , the variance of Hi vanishes as 
[H 2 ] - [H\ 2 oc 1/L d , with [HP] = (l/K)J2? =1 Hf. It 
is important that the variance decays with exponent d, 
i.e. there should be no critical exponent involved. The 
field Hi that maximizes T is just chosen to cancel the 
random field excess A# of Eq. I|14p. the square of which 
scales inversely with the volume. The behavior of the 
variance above T c is less relevant because here we no 
longer have phase coexistence, and so the OPD tends to 
a single Gaussian peak as L — > oo. The path lr, as well as 
the "equal-weight" path, then become meaningless any- 
way. 

IV. NUMERICAL TESTS FOR THE RFIM 

We consider the RFIM Hamiltonian of Eq. (|32p , using 
Gaussian random fields with a = 1.4, for which New- 
man and Barkema (NB) report as critical temperature 
T C NB w 3.6 [HI. Since the distribution of random fields 
is symmetric about zero, it also holds that H c = 0. The 
implications of modified hyperscaling will now be veri- 
fied. 



A. FSS using the symmetry path Is 

We first use the symmetry path. That is, we mea- 
sure PL,i{m) at fixed H = 0. Even though the criti- 
cal field H c = 0, spin reversal symmetry is broken in 
single samples, and so we do not expect PL,i(rn) to 
be symmetric (only after the disorder average [•] has 
been taken is the symmetry restored). To verify that 
symmetric distributions are rare, we consider the ratio 
A = J , Pj, ) j(m) dm/ J_ 1 Pi ! ,(m) dm. For a perfectly 
symmetric distribution A — 1/2 (the reverse is not nec- 
essarily true). Fig. [5ja) shows histograms of observed A 
values, measured at T = T C NB , and for various L. For 
each system size, K ~ 10, 000 random field samples were 
used; more details regarding this choice are provided in 
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FIG. 5. Investigation of the "typical" shape of PL,i(m) in 
the RFIM at T = T C NB using the symmetry path l s . (a) His- 
tograms H(A) for various system sizes L. Note the logarith- 
mic vertical scale, (b) Typical distribution Wl,% = In PL,i{m) 
for L — 14; for this distribution A ~ 1. 



the Appendix. It is clear that symmetric distributions 
are rare. Most distributions yield a value of A close to 
zero or unity, meaning that the "weight" is entirely con- 
centrated left or right of m = 0. Fig. [5jb) shows the 
logarithm Wl,i = PL,i(m) of one such "typical" distri- 
bution. A bimodal structure is revealed, but the peak 
heights are very different. If one plots Pi,i(m) itself it 
is clear that only a single peak survives. We conclude: 
by using the symmetry path, Pi, t i(m) is mostly a single 
peak. However, note that H(A = 1/2) is not zero: distri- 
butions whose "weight" is spread symmetrically around 
m = do occasionally occur. We return to this point 
later. 

The symmetry path Is does not lend itself well to 
measure free energy barriers, most distributions being 
single-peaked, but we can still probe sample-to-sample 
fluctuation^]. For each random field sample i, we calcu- 



A "work-around" to extract the barrier using the symmetry path 
can still be defined, see Appendix l"4l 



FIG. 6. Cumulant analysis of the RFIM using the symmetry 
path Is- For clarity, the curves are shifted by unity, and a 
logarithmic vertical scale is used, (a) (7i, con versus T for vari- 
ous L; the arrow marks T^ B . (b) The disconnected cumulant 
cA.dis versus T for various L. Note the absence of intersec- 
tions in {/i lC on, and their presence in t/i,dis- This is conform 
the modified hyperscaling scenario. 



late the magnetization (|m|), = \m\PL,i(Tn) dm, and 

the second moment {m 2 )i = m 2 PL i(m) dm, which 

arc then averaged to obtain [(|m|)] = (1/K) ^2^ = i(\m\)i, 
and so forth. In Fig. EJa), we plot the connected cu- 
mulant C/i 7Con (Eq. ([55)) ) versus T for various system 
sizes, while (b) shows the disconnected cumulant cTi dis, 
Eq. (|37[) . The striking result is that LTi^is reveals an in- 
tersection point, while L^rfim does not: exactly what is 
predicted by modified hyperscaling! From the intersec- 
tions in E/i )C iis) we conclude that the critical temperature 
is somewhat below T C NB . Above T c , the connected cumu- 
lant Ui )Con — > ir/2 as L — > oo. If one plots Ui tCon versus 
T for two values of L, an intersection will also be found, 
at some value Tl > T c ; sec for instance the curves for 
L = 14 and L = 16 in Fig.[6ja). As L increases, Tl will 
shift toward T c , but there is no intersection of U\ con at 

We still discuss the quenched-averaged distribution 
Qhim) = {l/K)^2f =1 PL.i(m), i.e. the arithmetic mean 
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FIG. 7. (a) The quenched- averaged distribution QL{m) of 
the RFIM obtained at T = T C NB , L = 14, and using the 
symmetry path Is- The peak-to-peak distance is proportional 
to the order parameter M, while the peak widths W reflect 
the sum of thermal and sample-to-sample fluctuations, (b) 
The leading cumulant of Ql(?ti) versus T for various L. 



of the individual (normalized) OPDs. Since PL i i(m) 
is mostly a single peak, located with equal probabil- 
ity at positive or negative values, Ql(to) is bimodal 
and symmetric about m = (Fig. E£a)). The pcak-to- 
peak distance corresponds to (twice) the order param- 
eter M = [(\m\)], but care is needed to interpret the 
peak widths W. The moments of Ql{jti) are of the form 
[{m k )\, and so the peak widths correspond to 

W 2 = [(m 2 )] - [(H)] 2 = Xcon /L d + Xdis /L d , (41) 

which is the sum of thermal fluctuations (set by Xcon) 
and sample-to-sample fluctuations (set by Xdis)- Conse- 
quently, the leading cumulant of Qi(m) becomes 



Ul,Q 



[(m 2 )] 



Xc 



L d M 2 ' 



(42) 



Using now the FSS expressions M oc L , Xcon oc L 1 ^ 
and modified hypcrscaling, we obtain 



Ui, Q - t/i, dis « L-^' 11 (T = T c 



(43) 



FIG. 8. (a) Variance [H 2 ]-[H] 2 of the "tuned" external fields 
Hi versus T using the path It for several system sizes, (b) 
The variance multiplied by L d . 



Hence, in the thermodynamic limit, U± q becomes iden- 
tical to C/i^dis- Plotting Ui : q versus T for different L one 
therefore also observes intersections (Fig. [7^1))). Note 
that, due to the correction term induced by the connected 
susceptibility in Eq. (T4"2"j) . one expects that for small L 
the intersections are more scattered than those for c/^dis; 
the data in Fig. |U[b) and Fig. [TJ^b) are compatible with 
this expectation. 



B. FSS using the path l r 

We now use the path lr, where for each random field 
sample i the external field Hi is tuned according to 
Eq. (T4TJI) . We first verify, in Fig. [51 that the variance 
of Hi indeed decays cx 1/L d . The raw data are shown in 
(a), while (b) shows the same data multiplied by L d . In 
the latter representation, the L-dcpcndcncc should can- 
cel for T < T c . This is confirmed by the collapse of the 
data of the larger systems; only the L = 8 data is some- 
what off, which indicates that this system may be too 
small for an accurate FSS analysis. The "swaying-out" 
of the curves at high T is a sign of entering the one-phase 
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FIG. 9. Investigation of the "typical" shape of PL,i(m) in the 
RFIM using the path lr- Shown are histograms H(Ui) for 
various L obtained at (a) T = T C NB and (b) T = 3.35. The 
histograms peak at Ui = 1 implying that most distributions 
are bimodal. Note the logarithmic vertical scales, (c) Typical 
distribution Wla — In Pt,i(m) obtained for L = 14 and T — 
T C NB . Since Wl,i is bimodal, a free energy barrier AFL,i can 
be extracted (vertical arrow). 



region, where the path lr becomes ill-defined. 

By using Zr , we expect that most distributions become 
bimodal for T < T c . In Fig.[9j we show histograms of ob- 
served cumulant values, for T = T C NB (a) and T = 3.35 
(b). We believe the latter temperature is closer to the 
true T c , based on the intersections of the disconnected cu- 
mulant. Fig. IDJb). The histograms peak at U\ = 1, con- 
firming that bimodal distributions dominate. An exam- 
ple distribution is shown in Fig.^c), from which a barrier 
AFls can be accurately extracted (vertical arrow). We 
remind the reader that the barrier is to be obtained from 
the logarithm of PL,i(m). Note also an important finite 
size effect in the histograms H(U{). At T = T C NB , for 
increasing L, a shoulder develops at Ui ~ 7r/2, mean- 
ing that single-peaked distributions become more likely 
in larger systems. In contrast, H(U{) at T = 3.35 reveals 
no such effect. We believe this indicates that T = T" C NB 
is actually above the critical temperature. Above T c , 
in the thermodynamic limit, the OPD is single-peaked. 
Hence, Fig. H^a) shows the evolution toward this shape. 
The convergence with L is clearly very slow, and much 
larger systems are required before single-peaked distribu- 
tions would dominate bimodal ones in finite-size simula- 
tion data. 

The path lr facilitates a first test of the scaling of the 
quenched-averaged barrier AFl = (l/-?0 A-Fl.i, 
where the sum is over all K considered random field sam- 
ples. For distributions where a barrier cannot be mean- 
ingfully defined, such as single or triple peaks, AFl^ is 
set to zero. In Fig.fTUTa), we show AFl versus T for var- 
ious L. Following modified hypcrscaling, we expect AFl 
to scale conform Eq. (|3"5|) . Hence, plotting L^ 6 AFl ver- 
sus tV-l", i = T/T c -l, the curves for different L should 
collapse, provided suitable values of 8, v, and T c are used. 
This result is shown in Fig. fTUT b). Here, = 1.5 was as- 
sumed, and by varying v and T c , a data collapse is indeed 
obtained (the plot uses v = 1.9 and T c = 3.32). We have 
verified that by using 9 = 0, i.e. the value of the pure 
model, no data collapse is obtained. The estimate of v 
is rather large, but still within the range of values re- 
ported in Table U Note also that T c used in Fig. ITUTb) 
agrees with that of the disconnected cumulant intersec- 
tions, Fig. [SJb). 

We now propose one additional method to locate the 
critical temperature. To this end, recall the FSS expres- 
sions AFl oc L e and Xcon oc L 1 ^ . Since 6 = ~f/v, the ra- 
tio K = Xcon/ AFl becomes L-independent at criticality. 
One can thus locate T c by plotting k versus T for vari- 
ous system sizes, and look for intersection points. This 
approach has the advantage that the critical exponents 
themselves need not be provided. The connected sus- 
ceptibility is obtained from the individual distributions 
P L ,i(m) using xcon = (L d /K)Y J f =1 {{m 2 ) l - (|m|)f). In 
Fig. fTTTa). we plot k versus T for various L. The data 
indeed intersect, providing important confirmation that 
the barrier scales with the same exponent as the con- 
nected susceptibility at criticality. For completeness, we 
show in Fig. Illf b) the disconnected cumulant C/^dis ver- 
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FIG. 10. FSS of the free energy barrier AF L in the RFIM 
using the path lr- (a) AFl versus T for various system sizes, 
(b) The same data as in (a) but using scaled variables. The 
validity of the scaling form, Eq. (|39p . is confirmed by the 
collapse of the data from the various system sizes onto a single 
curve. 



FIG. 11. FSS in the RFIM using the path l r . (a) k versus T 
for various L. Note the intersection point, which yields an es- 
timate of T c . (b) The disconnected cumulant, Eq. (|37[) , versus 
T for various L. The intersections again yield T c . For clarity, 
the cumulant curves are shifted by unity, and a logarithmic 
vertical scale is used. 



sus T for various L (now obtained using the path It). 
The curves also intersect, and do so remarkably close to 
the intersections of k. Based on Fig. QTJ we (VFB) re- 
port T C VFB w 3.315 ± 0.050, where the error reflects the 
scatter in the various intersection points (here the data 
of the smallest system L = 8 was ignored). 

We now turn to the distribution Vl{{\iti\)), defined 
as the probability of a particular random field sample 
yielding a magnetization (\m\). At criticality, we antic- 
ipate scaling of this distribution, conform Eq. (|36)l . We 
have explicitly measured T^x, ((|»n|)) by accumulating a 
histogram of (\m\) values at T = T C VFB using the path 
It- The resulting distributions are shown in Fig. [FZ[ The 
salient features are a sharp peak, and a long tail extend- 
ing to lower values. The fact that VL{{\m\)) features a 
sharp peak is consistent with C/^dis being close to unity 
at criticality. Since /? ~ in the RFIM, the scaling vari- 
able x = L~P'"(\m\) is identical to (\m\) itself, and so the 
"raw" distributions for different L should already over- 
lap with each other. Within numerical precision this is 
confirmed, but it is clear that the data in Fig. do not 




FIG. 12. r L ((\m\)) for the RFIM at T — T C VFB , using the 
path lr, and for various system sizes. 
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allow for any meaningful estimate of (3 /v. 

The point that we wish to make, however, is a differ- 
ent one. The fact that Pl((I to D) features a long tail 
means that occasionally a distribution -Pl.^to) is ob- 
served with a significantly lower magnetization. Since the 
scaling form, Eq. (j3"rJ)) . implies that VL{{\m\)) retains its 
shape irrespective of L, the fraction of these distributions 
does not vanish in the thermodynamic limit. It is con- 
ceivable that distributions from the "tail" of 'Pl((\^\)) 
are also shaped differently. For instance, consider again 
the histogram H(U\) at T c (Fig. EJb)). The histograms 
peak at U% = 1, so most distributions -Pl^to) are bi- 
modal. However, H(U\) also features a tail, so distribu- 
tions with profoundly different shapes, although rare, do 
occur. In particular, the tail in H(U\) allows for three- 
peaked distributions to be present (for which U\ = 3/2). 
Indeed, such distributions are observed, and have been 
interpreted to signify first-order transitions [46|, or new 
phases [47|. Our point is that the long tail of Vl ((|rn|)) 
and its scale invariance at T c (implied by modified hyper- 
scaling) also allows for the presence of three-peaked dis- 
tributions (without having to assume a first-order tran- 
sition, or the emergence of a new phase). 



V. WIDOM-ROWLINSON MODEL WITH 
QUENCHED OBSTACLES 

It was argued by de Gennes that a binary mixture un- 
dergoing phase separation inside a random network of 
quenched obstacles belongs to the universality class of 
the RFIM [48| . The argument is expected to hold when 
the obstacles display a preferred affinity to one of the 
phases. In case there is no such preference, the argu- 
ment does not apply [H, [50[. Previous simulations @ 
have already produced evidence in favor of de Gennes' 
argument. To provide further confirmation, in particular 
to test the scaling of the free energy barrier (Eq. ©), 
we consider in this section the Widom-Rowlinson binary 
mixture (WRM) [Hi]]. The model consists of unit diam- 
eter spheres, species A or B, which may overlap freely 
except for a hard-core repulsion between A and B par- 
ticles. The model is investigated in the grand-canonical 
ensemble, where the relevant thermodynamic parameters 
are the fugacities, za and zb, of the respective species. 

At high fugacities, the WRM can be in two phases: a 
phase rich in A particles (the A- phase) when za > zb, 
and a phase rich in B particles (the B-phase) when 
zb > za- Due to the model's symmetry under the ex- 
change of A and B particles, the phase transition oc- 
curs at za = zb- Hence, in line with the Ising model, 
a symmetry path Is for the WRM can be defined as 
Is '■ za = zb- The transition line ends in an Ising critical 
point, at fugacity za = zb = z cr i t , below which mixed 
states appear 52 , [53| . Note that the phase transition in 
the WRM can also be considered a liquid-gas transition. 
By integrating out the B particles, the WRM maps onto 
a single component fluid, interacting via a short-ranged 




FIG. 13. Cumulant plot for the WRM without quenched dis- 
order. Plotted is Ui versus zb for different L. The intersection 
yields z CI .it = 0.9377(4) and Ut = 1.228(5). 
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FIG. 14. Extrapolation of zl(x), x £ {x>X4 iXt}> according 
to Eq. (|45|) for the WRM with quenched disorder. From the 
linear fits z cr it = 0.9376(5) is obtained. 



attractive potential [5l|. The fugacity zb then plays the 
role of inverse temperature, the A-phase corresponds to 
the liquid (characterized by a high particle density) , and 
the -B-phase to a gas (low particle density). 



A. pure mixture 

We first consider the pure WRM, i.e. without quenched 
obstacles. We simulate using cubic boxes with periodic 
boundary conditions (see Appendix [5J . The analogue of 
the Ising model OPD is the distribution Pl{pa), defined 
as the probability for a system of lateral extension L to 
contain Na = PAL d particles of species A. Since we are 
ultimately interested in locating the critical point, only 
OPDs lying on the symmetry path lg are considered in 
this section, which leaves as the single free paramc- 
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tcr. Note that we could also have defined the OPD as 
Pl{pa ~ Pb)i thereby directly exploiting the symmetry 
of the WRM. However, most fluids lack such an obvi- 
ous symmetry, and by using Pl(pa) we ensure that our 
analysis remains generally applicable. 

Above the critical fugacity, zb > z CT it, Pl(pa) is bi- 
modal: the peak at low (high) density corresponds to the 
gas (liquid) phase. When zb < Zcrit> the OPD features 
a single peak, corresponding to a mixed state. The ana- 
logue of the magnetization is defined as m = pa — (pa), 
which is readily substituted in Eq. (fTS| to yield the order 
parameter and susceptibility. Additionally, we define a 
"generalized" susceptibility pjol l45j 

X4 ee L 3d «m 4 ) - 4(M)<H 3 ) + 12<m 2 )(M> 2 

-3<7n 2 ) 2 -6<|m|> 4 ). 1 ' 

The most straightforward method to locate the critical 
point is from intersections of the Binder cumulant for 
different L. For the pure WRM, we find that a sharp 
intersection of U\ can be found easily (Fig. [T3| . Another 
method to locate the critical fugacity is via the extrapo- 
lation of the finite-size extrema of x an d X4- In a finite 
system of size L, x reaches a maximum at fugacity zl(x)j 
which is shifted from z cr ;t as [45j 

^crit - z L ( X ) cx L- X l\ (45) 

with v the correlation length critical exponent. In addi- 
tion, X4 reaches a minimum and maximum, at respective 
fugacities z_l(x4 ) an d 2£,(x7), which are also shifted ac- 
cording to Eq. (|45|) . Hence, plotting Zl(x), z l{xI )> an d 
z^ixt) versus L" 1 /^, and then linearly extrapolating to 
L — > oo, three additional estimates of z cr jt are obtained. 
For this extrapolation, hyperscaling is not required, but v 
needs to be provided. In principle, v can also be taken as 
a fit parameter, but this requires data of extremely high 
quality. For the pure WRM, which belongs to the Ising 
universality class, v is known (cf. Table|I|. In Fig.QTl the 
extrapolation is demonstrated; the resulting estimates of 
z cr j t are similar and agree with the cumulant intersec- 
tions. Combining all results, we obtain z cr j t = 0.9377(5), 
where the error reflects the scatter between the individual 
estimates. This value is in good agreement with previous 
results [HHi. 
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FIG. 15. Cumulant plot for the WRM with quenched obsta- 
cles. Plotted is [/i, C on (Eq. (|35}) versus zb for different L. 
The intersection of system sizes L < 10 is attributed to a 
crossover effect from Ising universality to RFIM universality. 
Curves for L > 10 no longer intersect at this point indicating 
that for these system sizes the crossover to RFIM universality 
has largely completed. 
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FIG. 16. Variation of k with zb for various L for the WRM 
with quenched obstacles (analogue of Fig. Illf a)). The in- 
tersection of the curves indicates a critical point around 
zb ~ 1.35. 



B. mixture with quenched obstacles 

We now consider the WRM model with quenched ob- 
stacles. We use spherical obstacles, species X and Y, 
having the same diameter as the (mobile) A and B par- 
ticles. The total number of X and Y obstacles equals 
Nx = Ny = Pq x L d , rounded up or down at random to 
the next integer. The obstacles are distributed randomly 
at the start of the simulation, irrespective of overlap, af- 
ter which they remain quenched: this defines one disorder 
realization i. Next, A and B particles are introduced, and 
grand canonical MC is used to construct Pl,i{pa) for that 



disorder realization (see Appendix) . The ^4-particles (B- 
particles) have a hard-core interaction with A-obstacles 
(Y-particles) but do not interact with F-obstacles (X- 
obstacles) . The original motivation for this choice was to 
restore the symmetry line Is ■ za = zb in the disorder 
average. However, in what follows, we will use the path 
It, whereby za is tuned for each realization of disorder 
such that d(pA)/d\ogZA is maximized. 

We still need to specify the obstacle concentration pq. 
For a noticeable random field effect, the thermal correla- 
tion length £ should be large compared to the typical dis- 
tance £q between obstacles. Following the FSS "Ansatz" 
£ cx L, this implies L £q. If L is too small, crossover 
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FIG. 17. Scaling plot of the free energy barrier, according to 
Eq. (|39[) . for the WRM with quenched obstacles (analogue of 

Fig. nab)). 




FIG. 18. Extrapolation of zl(x), x £ {[x], [xi]~ j [X4] + }, via 
Eq. (|45|) for the WRM with quenched disorder. The plot uses 
v — 2.1 (as obtained in Fig. I17[). and predicts z cr it = 1.37(2) 
in the limit L — > oo. 



scaling is observed (in this case from pure Ising to ran- 
dom field Ising Q). From these considerations, choosing 
a high value of pq seems optimal. The disadvantage is 
that also z cr j t will then be very high, which makes the 
grand canonical MC approach inefficient due to a high 
particle density. Clearly, a compromise needs to be made: 
we use pq = 0.02. This value is small compared to typ- 
ical density of the mobile species, e.g. pa = Pb ~ 0.38 
at criticality in the pure WRM, and certainly is below 
the percolation threshold; we thus remain in the limit 
of weak random fields. For the chosen obstacle density, 
crossover effects are still strong in small systems. This 
can be inferred from Fig. [Jj)J where the connected cu- 
mulant Ui tCon versus zb for various L is plotted. The 
curves for L < 10 reveal an intersection point: this would 
be consistent with a conventional critical point featuring 
standard hypcrscaling. However, for L > 10, the inter- 



section has vanished, indicating that the crossover has 
largely completed. In what follows, we therefore discard 
the data for L < 10 in some of the analysis. 

Investigations involving disconnected quantities re- 
quire enormous simulational effort to generate data of 
sufficient quality (see Fig. [T5] in Appendix P2J) ■ For the 
WRM, an analysis of f/^dis along the lines of Fig. ITTT b) 
was not feasible. We therefore focus on the free energy 
barrier. We evaluate the distributions Pl,i{pa) along the 
path It , and for each distribution, we "read-off" the bar- 
rier, which is then averaged over the samples to obtain 
AFl. Wc first consider the variation of k versus Zb for 
different L, i.e. the analogue of Fig. [TlTa). This data 
is shown in Fig. 1161 from the intersection we conclude 
that the critical fugacity is around zb ~ 1.35. To get 
the critical exponents, we consider the scaling of the free 
energy barrier. Assuming RFIM universality, the varia- 
tion of AFl with zb should follow Eq. (|31?|) . where now 
t = (z cr it — zb)/ 1 zb- In the vicinity of zb ~ 1.35, i.e. as 
indicated by Fig. [16j we indeed find that a collapse of 
the curves can be realized for = 1.32, v = 2.1, and 
Zcrit = 1-37 (Fig.QTJ). As a consistency check, we attempt 
to obtain z cr j t from the extrapolation of the extrema of 
the susceptibilities using Eq. (|45j) . The observables x an d 
Xi of the pure model are now replaced by their disorder- 
averaged counterparts [x] and [xi], and v = 2.1, i.e. the 
estimate from Fig. II Tl is used. The extrapolation works 
reasonably well (Fig. [THj) and for the critical fugacity we 
obtain the same estimate as before: z cr jt = 1.37(2). 



VI. SUMMARY 

Modified hypcrscaling, Eq. @, which is believed to 
describe systems belonging to the universality class of 
the RFIM, gives rise to rather unusual finite size effects 
at critical points: neither the order parameter distribu- 
tion, nor the free energy barrier AFl of interface forma- 
tion, are scale invariant. As a result, "standard" tech- 
niques to locate critical points, such as the "cumulant 
intersection method" [l4[, or the Lee-Kosterlitz method 
[43^ . break down. However, by carefully considering the 
consequences of Eq. ([2]) , alternative techniques to derive 
T c in random field systems can be derived. In this pa- 
per, we have proposed two such techniques. The first is 
based on the order parameter fluctuations between dis- 
order samples: modified hypcrscaling predicts that these 
are scale invariant at T c . This property can be used to 
locate T c by measuring the disconnected cumulant t/i.dis 
(Eq. (|37p ) versus temperature for various system sizes: 
at T c , curves for different L intersect. Indeed, simu- 
lation data of the RFIM confirm the scaling of Ui^is 
(Fig. Eljb) and Fig. ITTTb)). In contrast to conventional 
critical points, there is no intersection of the connected 
cumulant Ui tCon (Eq. (|35]l) in the RFIM at T c . However, 
in small systems, there may be crossover effects. In this 
case, an apparent intersection in Ui. con is observed, at 
T > T c , but it vanishes in larger systems; such was the 
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case for the WRM (Fig. ITS]). 

The practical disadvantage of measuring t/i.dis is that 
many disorder samples must be averaged over if mean- 
ingful results are to be obtained. Particularly for more 
complex systems, such as off-lattice fluids, an economic 
alternative is to consider the free energy barrier AFr, of 
interface formation. Due to modified hyperscaling, the 
barrier diverges AFr, °c L e at T c , with 9 the violation of 
hyperscaling exponent. The consequences of this diver- 
gence are easily detected in simulations, as was demon- 
strated for the RFIM (Fig. [TO] and Fig. [TQX)), and the 
WRM (Fig. [H] and Fig. H7J) . In case of the RFIM, the 
estimate of T c obtained from the scaling of the barrier 
was fully consistent with that obtained from the inter- 
sections of U 1: dm (Fig. E]). Our results for the WRM 
provide further confirmation that fluids with quenched 
disorder indeed belong to the universality class of the 
RFIM, consistent with the conjecture of de Gennes [Hf. 

We have also commented on the variations in shape of 
the order parameter distribution between samples. There 
is some question as to whether distributions with three 
peaks signify first-order transitions fit! ] , or the emergence 
of new phases 47 1 . Our view is that modified hyperscal- 
ing also allows for these shape variations. While our data 
indicate that at T c , and using the path lr, the majority of 
distributions is bimodal, a fraction of distributions with 
different shape is not ruled out (Fig. M^>) and Fig. DPT 



Finally, we remind the reader that the divergence of 
the free energy barrier at T c will also influence the dy- 
namics. Taking the RFIM with single spin-flip dynamics 
as example, it follows that the largest relaxation time in 
a finite system at criticality is given by an Arrhenius' 
type formula 



lnrocZ/ (T = T C ) 



(46) 



This is in contrast to the pure model, where the relax- 
ation time (not its logarithm) scales r oc L z , with z the 
"dynamical critical exponent" . Such a power law for the 
logarithm of the relaxation time is the hallmark of "ac- 
tivated critical dynamics" . In fact, if we are somewhat 
above T c , but the system size is still less than the corre- 
lation length, L < I, Eq. (05]) still holds! As L > £, the 
system size in Eq. (|46|) gets replaced by £, and we recover 
Eq. (jj), as proposed by Villain [2f| and Fisher [H]. A 
direct study of the dynamics of a kinetic version of the 
RFIM would be illuminating, but goes beyond the scope 
of the present paper. 
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Appendix: Simulation Details 



1. Wang-Landau sampling 

The simulations of the Ising and RFIM were performed 
using Wang-Landau (WL) sampling [H-jj . The OPD is 
written as 

PlAS) oc g L ,i(S)e HS/kBT , (A.l) 

with S = L d m the total instantaneous magnetization, 
and gL,i(S) some generalized density of states (DOS). 
Note that the DOS depends on system size L, tempera- 
ture T, and random field sample i, but not on the external 
field H (for the pure Ising model, there is no dependence 



on i either, of course). At the start of each simulation, 
we generate a sample of random fields i. We then per- 
form single spin-flips, whereby one of the spins is chosen 
at random, and its orientation reversed. Let the total 
magnetization and energy at the start of each spin-flip 
be given by So and Eq, respectively, and afterward by Si 
and E\] each spin- flip is then accepted with probability 

a{S ,E ^ S u Ei) = 

Note that the energy above refers to the configurational 
part of the Hamiltonian only, i.e. the nearest-neighbor 
interaction and the coupling to the random field, but not 
the coupling to the external field. 

The DOS is a-priori unknown, and is initially set to 
unity gL,i{S) = 1. After each attempted spin-flip, one 
"updates" the DOS g L ,i( s ) ~> f x 9lA s )i witn s thc 
magnetization of the system after the attempted spin- 
flip. The update is performed irrespective of whether the 
spin-flip was accepted; the initial modification factor / = 
e w 2.72. We also update a histogram h(S) — > h(S) + 1, 
counting how often a state with magnetization S was vis- 
ited. This procedure is repeated until h(S) has become 
sufficiently flat, which completes one WL iteration. We 
use the criterion (h max - h min ) / '(/i max + h min ) < 10~ 5 , 
with h m in and h max the smallest and largest entries in 
h(S), respectively. After the first WL iteration, the mod- 
ification factor is reduced / — > f 1 ^ 2 , the histogram h(S) 
is reset to zero, and the procedure is repeated. WL iter- 
ations are continued until / has become small such that 
changes to the DOS become negligible. For each DOS, 
we typically performed 150-250 WL iterations. Once the 
DOS is known, the OPD can be calculated for arbitrary 
values of H using Eq. (|A.1|) . 

2. Importance of disorder averaging 

To accurately determine disorder averages, the OPD 
PL,i{ m ) ' s measured i = 1,...,K times. In particular 
disconnected quantities require a large number of disor- 
der samples if meaningful results in the critical regime 
are to be obtained. Fig. Q1J] shows a typical "running av- 
erage" of Xcon and Xdis versus K. While Xcon saturates 
to a plateau already after 1000 samples, the convergence 
of Xdis is noticeably slower. The data of Fig. indicate 
that K should exceed several thousands at least. Away 
from the critical point, Xdis is no longer divergent, and 
here we expect that lower values of K will also suffice. 

3. Histogram reweighting in temperature 

A key ingredient in this work is the use of histogram 
reweighting in the temperature-like variable. We perform 
our simulations at only a few distinct temperatures, and 
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10000 



FIG. 19. "Running averages" of the connected and the dis- 
connected susceptibility versus the number of random field 
samples K. The data were obtained for the RFIM using the 



symmetry path Is, L = 
rithmic horizontal scale. 



14, and T = IT 



Note the loga- 



extrapolatc to other values using histogram reweighting 
[55[ | . This requires that the joint two-dimensional prob- 
ability distribution Pl,i(S,E), of the magnetization S 
and energy E, is known. Again, as in Eq. (|A.2[) . E refers 
to the configurational part of the Hamiltonian only. If 
PL,i(S,E) is measured for T = To and H = Hq, it can 
be extrapolated to other values using a generalization of 
Eq. CO} 



PL,i(S,E)\ TuHl <xP L<i (S,E)\ TotHo x e 



ShS-SpE 



(A.3) 



with Sh = (Hi - H )/k B Ti and 5f3 = l/k B Ti - 
1/fcsTo. The practical problem is that two-dimensional 
histograms require considerable disk space, which in the 
case of quenched disorder is multiplied by a factor K. 
Fortunately, an excellent approximation can be used to 
drastically reduce storage requirements @ . Without loss 
of generality, we write 



P L 4S,E)=P L 4S)xgf}(E), 



(A.4) 



where gf\(E) is the probability distribution of the en- 
ergy measured at states with the same magnetization S. 
The approximation is to assume that g L ■ (E) is Gaussian, 
and so is fully specified by its first two moments. For each 
random field sample, the two-dimensional histogram of 
Eq. (|A.4[) then requires only PL,i(S) to be stored, plus 
the "functions" {E) L ,i(S) and (E 2 ) La (S). 



4. Alternative method to measure the barrier 

It is also possible to measure the quenched-averaged 
free energy barrier AF^ using the same external field for 
all samples. To be concrete, consider the OPD PlA^) 
of the RFIM obtained at fixed H = 0, i.e. the symmetry 




FIG. 20. Demonstration of an alternative method to extract 
the free energy barrier AFl. The data in this figure were 
obtained for the RFIM using L = 14. (a) The quenched- 
averaged free energy distribution Wz,(m) constructed with 



the recursion relation of Eq. 



m at T = T c v 



a free en- 



ergy barrier AFl can be meaningfully extracted, (b) The 
corresponding variation of AFl versus T, compared to the 
"original" method, where AFl is averaged over individual 
samples using the path lr- 



path, with total magnetization S = ~L d , —L d + 2, . . . ,L d 
and i = 1, . . . , K. We define the quenched-averaged free 
energy difference between "adjacent" states as 



AW L (S-2,S) 



1 K 
K ^ 



PlAs-z) 



(A.5) 



which can be used to construct a total free energy Wl (S) 
by means of recursion 



W L {-L d ) =0, 

W L {S) =W L {S-2) + aw as - 2, S). 



(A.6) 



Fig. l2Tira) shows the typical shape of the free energy ob- 
tained in this way for the RFIM. The distribution is bi- 
modal, and a free energy barrier AF^ can be meaning- 
fully "read-off". As it turns out, this barrier is very simi- 
lar to that obtained by averaging over individual samples, 
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i.e. as was done in Fig. HOT a) using the path lr; a compar- 
ison is provided in Fig. I2"fjrb'). In fact, if one uses Wl(S) 
to perform the scaling analysis of Fig. [TCT b). excellent 
data collapses are also realized. 

An analysis in terms of Wl(S) is numerically conve- 
nient because extrapolations in the field variable H can 
be performed after the quenched average has been taken 

W L (S)\ Hl = W L (S)\ Ho + (F x - H Q )S/k B T. (A.7) 

This is particularly useful for fluids, where the critical 
field (chemical potential) is generally not known before- 
hand. However, it is not obvious what the peak positions 
and widths in Wl(S) correspond to. Based on our previ- 
ous work cumulants of e Wh ^ do intersect at T c , but 
(in hindsight) we believe it is safer to perform the cumu- 
lant analysis using the individual OPDs (as was done in 
this work). 



5. Simulating the Widom-Rowlinson model 



of species x is placed at a random location; the resulting 
new state is accepted with probability 



a(N x 



N x 



1) = min 



1. 



z x V 



1 



h 



(A.9) 



States with hard-core overlaps and states where Na is 
outside the window bounds are always rejected, irrespec- 
tive of the accept probabilities. 

While simulating in window Wk , we keep track of two 
counters, CZ and Ct ■ These count, respectively, how 
often the state with Na = k and Na = k + 1 was visited. 
From these counters, we construct the relative probabil- 
ity of these states via 



PL,,(fc + l) 



ri-Jk ■ 



(A.10) 



Having at hand this ratio for all windows Wk, the full 
distribution is constructed recursively 



We measure Pl,i{pa) using grand canonical MC and 
successive umbrella sampling (SUS) [56|. The simula- 
tions are performed in a periodic 3D cube of volume 
V = L 3 . At the start of each simulation, a disorder 
realization i is generated by distributing the obstacles X 
and Y randomly in the cube, i.e. the obstacles are allowed 
to overlap. We then perform grand canonical MC moves 
consisting of the insertion and removal of single A and 
B particles. In SUS, the full density range of interest is 
split into overlapping windows Wk- In the first window, 
Na is allowed to fluctuate between and 1, in the second 
window between 1 and 2, or, more generally, in the fc-th 
window Wk ■ Na G {k, k + 1}. There is no restriction on 
the number of B particles: Nb thus fluctuates freely in 
each window. 

For Na = the B-particles are an ideal gas in 
the volume allowed by the quenched Y"-particlcs so an 
initial state for Wq is easily constructed: we draw 
a number N from a Poissonian distribution P(N) — 
£ -zbV (zbV) n /N\, and randomly insert this number of 
J3-particles into the system, discarding all B-particles 
that overlap with ^-obstacles. As starting state for the 
subsequent windows Wk (k > 0), we take the last state 
of the window Wk-i preceding it, and equilibrate this 
state briefly for ~ 10 5 MC steps within the bounds of 
the new window. This works well in practice because the 
windows are small. 

The production run of each window Wk is performed 
using ~ 10 7 MC steps. Each step first selects a species, 
x G {A, B}, with equal probability. Then, with equal 
probability, the insertion or removal of a particle of 
species x is attempted. In case of removal, a particle 
of species x is picked at random and removed from the 
system; the resulting new state is accepted with proba- 
bility 



n a -i 



PlANa)<x ' 



a(N x — >• N x — 1) = min 

The factor fk is one when x = 
specified later. In case insertion is chosen, a new particle 



i ^r 1 

'z x V h 
B: for x = 



(A.8) 
A, it will be 



PL,i{k) 



N A -1 r+ 

II %r k \ (a.ii) 

fc=0 °fc 



where the proportionality constant follows from normal- 
ization. Note that Pl^Na) above is, of course, fully 
equivalent to the OPD Pl,i{Pa) that we wish to find. 

We now specify the factor fk for moves involving A- 
particles. Assuming a constant number of steps per win- 
dow, Eq. (|A.10[) suggests that optimal results are ob- 
tained when fk is chosen such that and C7 arc 
roughly equal, i.e. fk = PL.i(k)/ Ph,i{k + 1), which is the 
sought-for result itself. For the first window k = we use 
the pure model's optimal weight, which can be calculated 
analytically. For the subsequent windows, we linearly ex- 
trapolate Pla{Na) to calculate fk = Ph,i{k— l)/PL,i{k) 
to be used in that window. In practice, this choice is al- 
ready quite good, and the counts in the upper and lower 
bin consistently lie within 1% of each other. 

In view of the huge amount of disorder realizations re- 
quired, a mechanism that allows for histogram reweight- 
ing of results obtained at (za, zb) to nearby parameters 
(za, Zb) is indispensable. To facilitate this rewcighting, 
the joint probability distribution Pl,i{Na,Nb) is stored 
in compact form as described in Appendix QJ] the re- 
sults of that section trivially transfer to the WRM if one 
identifies S O N A and E <-> Nb- For the WRM with 
quenched disorder, the range m zb over which one can 
reliably extrapolate is too small to cover the full region of 
interest. In particular, simulation data obtained at the 
fugacity zl([x]) of the susceptibility maximum could not 
be extrapolated to the critical fugacity z ct it- We there- 
fore created two data sets per system size. One set with 
K = 2000 disorder realizations at zb ~ z l{[x\) used for 
locating the extrema of \ and Xi (Fig- HH) an d one set 
with K — 10000 realizations around zb ~ 1.4, which 
is close to Zcrit, for investigating the free energy barrier 
(Figs. Dl and [TJJ. 



